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Abstract. An iterative scheme is presented to solve analytically the relativistic fluid dynamics equations. 
The scheme is applied to longitudinal expansion, transversal symmetric and transversal asymmetric 
(triaxial) expansion as well. Within this scheme it is possible to describe the dynamics of a strongly 
coupled (i.e. conformal) medium for parameters referring to heavy- ion collisions at LHC. 
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1 Introduction 

■ Fluid dynamics became a standard tool for modeling the expansion dynamics of the strongly coupled 
quark-giuon plasma (sQGP) produced in ultrarelativistic heavy-ion collisions. For given initial conditions, 
e.g. inferred from transport models or Glauber type Monte Carlo simulations, the sQGP's evolution is 
followed in space and time until the density becomes so low that a description by hadronic degrees 

pi I of freedom is appropriate and the late history until and including the freeze-out can be described by a 

1^ , hadronic transport model. In fact, various observables, such as transverse momentum spectra, differential 

■ elliptic flow and higher flow modes as well, can be successfully explained for RHIC and LHC beam 
energies [1112. The initial conditions for the subsequent fluid dynamical evolution might be smooth, i.e. 
emerging from an average over many events, or might be fluctuating, i.e. corresponding to an event-by- 
event evolution [3] . The influence of initial state fluctuations on final state observables currently receives 
much attention (for a selection of recent works c/. [3]) which to a large extent is motivated by a detailed 
understanding of the particle distributions and correlations beyond the elliptic fiow mode as well as 

jy-^ ■ transport properties of the expanding medium. 
1 1 At the heart of the fluid dynamical description are the equation of state, accessible in lattice QCD 

■ calculations [3 [6l [71 18] , and transport coefficients, being less directly accessible by lattice QCD [9l ITOl 
! [H] due to their real-time nature. The AdS/CFT correspondence, however, allows estimates of the 

transport coefficients (cf. [11] [T^l [13] ) and further properties of the matter in the strong coupling limit, 
thus complementing perturbative weak-coupling approaches. In such a way, quantitative results for the 
properties of the sQGP can be obtained. Nevertheless it remains not yet clear, to which extent AdS/CFT 
results can be adopted. Intriguing questions concern the details of the confinement transition region. 



^ ■ position of the conjectured (critical) end point in the phase diagram and the temperature dependence of 

the transport coefficients in various regions of the phase diagram jl4| . 

Relativistic fluid dynamics might be grouped into perfect fluid dynamics, i.e. neglecting dissipative 
phenomena, dissipative fluid dynamics, e.g. as Landau-Lifschitz formulation [151 116] . and extensions 
thereof such as the Israel-Stewart type setup [T7] or transient fluid dynamics [IH] . These fluid dynamical 
descriptions {cf. [191 [20] for further discussion of relativistic dissipative hydrodynamics) represent partial 
differential equations which are usually solved numerically by employing suitable differencing schemes. 
Only in case of high symmetries of the flow, analytical solutions are at hand. Examples, relevant for 
heavy-ion collisions, are Bjorken's flow [^T], Gubser's flow [12] and the flow patterns studied by Csorgo 
et al. [221 [m [IS]- (Further flow patterns refer to Landau type initial conditions, cf. [2S]-) With respect 
to the important role of elliptic flow, however, more realistic solutions are desirable. The elliptic flow 
analysis with account of longitudinal pressure gradients calls for the time evolution of an initially triaxial 
matter distribution. This goal is envisaged in the present paper. We are going to utilize a scheme 
for solving algebraically the fluid dynamical equations for arbitrary but smooth and analytically given 
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initial conditions. We test an iterative procedure which ahows for solving fully analytically the fluid 
equations of motion up to a given order in the expansion of powers of the time coordinate. Such iterative 
procedures have been established in the context of general relativity (c/. [H]). In employing the AdS/CFT 
correspondence, the iterative procedure is exploited, e.g. in [28j, to solve the five-dimensional vacuum 
Einstein equations with negative cosmological constant and to derive useful insights in the behavior of 
the four-dimensional finite-temperature, strongly coupled medium in the conformal limit. 
Here, we outline such a solver for the four-dimensional fluid dynamics in section[2]and make it explicit for 
ideal fluid dynamics in section |31 Various applications up to the triaxial expansion dynamics including 
eccentricity as basis of the elliptic flow for conditions referring to heavy-ion collisions at LHC are presented 
in section |4l We discuss the limitations of the method in section [5] and draw conclusions in section |6l 



2 Description of the iterative scheme 

The basic idea of the iterative scheme to be described in the following paragraphs is to transform the 
differential equations of interest (in our case: local energy-momentum conservation) into an infinite set 
of algebraic relations between the Taylor coefficients of the unknown functions, which can be solved iter- 
atively. For achieving this, the unknown functions need to be formally Taylor expanded and afterwards 
the series is plugged into the differential equation. Comparing coefficients of the terms with the same 
power in the time-coordinate leads to the desired set of algebraic relations which is solved iteratively up 
to a truncation order. 

For this purpose, let us assume that the only fluid dynamical equations governing the dynamics of the 
sQGP are 

T^"'.^^ EE d^T^"" + r£:^r"'' + rfl^T^'^ = 0, (1) 

where the Christoffel symbols Fj^^ take into account non-Cartesian coordinates mapping out the four 
dimensional Minkowski space-time. In the absence of conserved charges the symmetric energy-momentum 
tensor T^'' is assumed to be related by a constitutive equation to energy density e, pressure p, four- 
velocity and additionally to gradients of these quantities in the form 

Tf"" = ew^w" + pA^"' + , (2) 

with A^"^ = u'^u'^ — g'^'^^heing the projector orthogonal to u'^. The quantity 11^'^ uncovers the dissipative 
effects. An equation of state e = e{p) and the normalization w'^m^ — 1 closes ([Ij to a system of four 
evolution equations for four independent quantities, e.g. e and the spatial components it" of the four- 
velocity. The Landau Lifschitz condition it^II'"^ = provides a unique definition of the velocity (c/. 
[ini for that issue). Cauchy initial conditions are specified by giving e{6o,S^) and u"(0o,O on the 
hypersurface 6 ^ 9o for a suitable time coordinate 9. For the subsequent calculations we use the 



coordinates ^'^ ^ (9,77,a;_L) being related to Cartesian coordinates {t,x) via 9 — In^yt'^ — (x^)"^ , 

= ^ In l^p- with tq being a reference proper-time; we choose it as the instant where the hydrodynamical 
evolution starts. The transverse coordinate vector x± = {x'^,x^) is equal to the Cartesian one. In these 
coordinates, the non-zero Christoffel symbols are 

r?! ^ = T^i = FJq = 1, (3) 

and the metric tensor g^^i, is c?iag(Tg e^®, — Tq e^®, — 1, — 1). The constitutive equations can be used to- 
gether with the equation of state and the normalization condition of the four-velocity to obtain the 
components of T^'^ as functions of four independent variables. We choose these variables to be T"", i.e. 

In principle, four arbitrary functions can be used to parametrize the components of the energy-momentum 
tensor. Since we want to apply the iterative scheme to relativistic hydrodynamics, one could choose the 
energy density and the spatial components of the four-velocity as well. There are two reasons why we 
not do so. First, setting up the iterative scheme is most straightforward for this choice and second, 
we want to keep the discussion in this section general, so the iterative scheme is not only applicable for 



^In Minkowski space we use the "mostly minus" sign convention for the metric. 
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ideal hydrodynamics, but also for other systems obeying local energy-momentum conservation. However, 
there are some subtleties concerning energy-momentum tensors which contain derivatives as in the case 
of viscous hydrodynamics. For such a situation Q needs to be modified (see appendixM . One important 
remark is that, when the setup is specified {e.g. by choosing to investigate ideal hydrodynamics with a 
certain equation of state), the functional form of /^^ is fixed and all its derivatives with respect to 
are given. Equation ([1]) can be rewritten as 

doT^^« = -daF'' ~ KpfP" - KJ^". (5) 

Latin letters for indices refer to spatial coordinates ranging from 1 to 3. Since /'^'^ depends solely via 
its arguments on the space-time coordinates, the spatial derivative of must be transformed into 
spatial derivatives of T"°. 

The iterative scheme we want to describe is based on the Taylor expansion 

oo 
fe=0 

where we have chosen without loss of generality coordinates for which the initial conditions are given at 
the manifold 9 = 0. Analogue expansions hold for the Christoffel symbols and for 



oo oo ^ 

k=0 ■ k=0 

Here f^j^.^ = '§gkf^'^iT"'^{0,^))\g^o- The time derivatives of Z'^" are then transformed according to the 
chain rule into time derivatives of T"'^ (see appendix [B|. Inserting these expansions into ([5]) yields 

oo ^ oo _j 

k=0 ■ ^ k=0 



oo /c 

2^2^ (k- l)W. V {k-l>phl) ^ {k-l)auJ(l) ) 



(k-iy.i 

It is important to note that /^'^■j contains only time derivatives of T"'^ up to order m. By comparing 
the coefficients, ^ can be decomposed into a series of algebraic recurrence formulas relating the Taylor 
coefficients: 



(=0 ^ ^ 



pp 

^ {k~l)vpJ(l) 



f(l) +^{k-l)a^f(l)) _ (9) 



Accordingly, the (fc -I- l)th Taylor coefficient T^i^_^^^ is algebraically to be calculated by lower-order 
coefficients. This makes an iterative solution possible. The initial conditions determine T^^. Using 
the equation with k = one can calculate from that. Applying the equation with k = 1 one can 
calculate T/^? etc. In such a way, arbitrarily high Taylor coefficients are calculated iteratively. 



= T"^{6 = 0,^) is supposed to be a smooth (differentiable) function. The determination of high- 



(2) 

-(0) 

order time derivatives is therefore via © reduced to calculating high-order spatial derivatives of T°'^{6 = 

For the (8, 77, xj_) coordinate system, the Christoffel symbols ^ are constants. Therefore only the terms 
with k = I contribute, and (O reduces to 



rpfJ.0 9 



QCa-lik) {Q)iypJik)^ ^ {0}auJ{k)) ■ ^-^"^ 
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3 Ideal fluid dynamics 



The scheme described in the preceding section is, in principle, apphcable to any physical system obeying 
local energy-momentum conservation. Therefore, it can be used to calculate solutions of the relativistic 
hydrodynamical equations. For ideal fluid dynamics, where XI^" in © is neglected, the relation ^ 
becomes explicitly 

with 



/ = 



gooT™' ) A (12) 



A = gabT^'^T''^ (13) 
(<7ooT"o' + a) goo 



By definition /^'^ — /"^ = T^". The four- velocity components are given by 



2gooTOO^ -A + T00^J-5A + AT^^^goo 



= ^ , (15) 

2^0^ A + Too'^oo 



with 



5A + 2.gooro°' - T^^Jg^o^-^A + 4gooTO" 



2VIV5ooroo' + A 



(17) 



The above equations hold for the conformal equation of state e = 3p, but can be generalized to others 
in a straightforward manner. They are obtained by considering the (/^O) components of These four 
equations are solved for the pressure p and the spatial components of the four velocity u° applying the 
equation of state and the normalization of the four velocity. Therefore, one can express p and as 
functions of the components of the energy-momentum tensor (c/. (ITil) and (HH)). The expressions 
p(r"°), w''(T"°) are plugged into ^ from which the algebraic relations (IT^ can be read off by applying 
(HI). We note that the equations (|TT|) - (fT7|) apply for any time-orthogonal coordinate system. 



4 Examples 

The iterative scheme is implemented in MAPLE to calculate derivatives and consecutively the expansion 
coefficients. 



4.1 Longitudinal expansion 

Let us consider first a purely longitudinal expansion. Following [30j we employ for the initial conditions 

2 

at 80 = a Gaussian energy distribution, e(0o, if) = eo exp{ — ^^-j with a = 3.8 for LHC energies and a 
synchronized flow y(9o,?7) — i]. The latter one is important for the flnding, already anticipated in |30| . 
that the energy density evolution can be approximated by e(0,r/) w e(0o, ?y) exp{— 10} and the flow 
velocity is only modified on the 2% level by a time dependent deviation from the Bjorken flow pattern 
characterized by a(0) according to y = 77(1 + a(9)). This statement is made more explicit below. 
Figure [T] exhibits e{Q,r])/e{<do,V — 0) (a), 6(6,77)/ e(9o,?7 = 0)exp{|6} (b) and y{Q,i]) — 77 (c) as a 
function of rj for various time instants. The iterative scheme has been used up to order N = 48. In figs. 
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[T](a) and (c), the iterative results (curves) are compared to a numerical solution of the hydrodynamical 
equations (c/. eqs. (10) and (11) in [5(7) obtained with a modified forward integration scheme (symbols). 
One observes agreement of both solution methods up to = 3.2, i.e. t = roe® = 24.5 tq. At later proper 
times r, the truncated iterative result shows some instability types, signaled by irregularities in e{r]) 
accompanied by negative energy densities at large values of rj. However, these can be separated from the 
physically interesting region. For larger values of Q, e.g. Q sa 4.0 {i.e. r w 55 tq), these instabilities even 
influence the energy density at the center. But at this proper time, the medium has cooled down so far 
that hadronisation is almost complete. In fig. [T] (c), the value of a{&) can be read off as the slope of 
the appropriate curve. Even for the largest depicted time Q = 3.2 the slope A(y — ?7)/A?7 being equal to 
a(6) is found to be small: 0.25/12 « 0.02. Since the slope of y{i]) increases monotonously with Q, the 
absolute value of a{<d) for < 6 < 3.2 is always smaller than the above mentioned 2%. The conclusion 
for this simplified flow pattern is that the Bjorken type scaling of the energy density oc exp{— 10} and 
synchronized fiow y — 77 « are excellent approximations for mid-rapidity \r]\ < 0.9, accessible in the 
ALICE experiment. Analytical corrections to the Bjorken behavior that can be obtained utilizing the 
iterative scheme are discussed below. 

Taking the initial temperature as T(6o,?7 ~ 0) ~ 430 MeV corresponds to the energy density e{Qo,r] « 
0) ^ 65GeVfm~^ (obtained with the s95p-PCE equation of state in for lattice data see e.g. [32]). 
Following the evolution until freeze-out at Tf.o. ~ 165 MeV, corresponding to the energy density of et.o. ~ 
500MeVfm~^, requires to go up to 6 = ln(T/To) = ln(e(6o)/et.o.)^^^ ~ 3.7. Figure [T] demonstrates that 
the iterative scheme, truncated at A'' = 48 is stable in this time interval. 

To expose the importance of the initially synchronized Bjorken type fiow pattern we chose y(Oo,77) — 
(1 + a{Qo))il- Foi' values of a(0o) not to far away from the Bjorken flow {i.e. a(6o) = — 1 . . .0.5), the 
deviation of the energy density from the reference system (a system with initially synchronized fiow, i.e. 
a(0o) = 0) can be parametrized by the series 



ercf 



^d„e"e-''"®. (18) 



Here, e denotes the energy density of the system with initially not synchronized flow, and Crof is the 
energy density of a reference system with the same initial energy distribution, but Bjorken type initial 
fiow. The parameters in this expansion can be directly linked to the Taylor coefficients of the energy 
density, by differentiating ([TS]) and evaluating the result at Qq. Apart from constant factors, the left 
hand side is the difference of the Taylor coefficients of e and Crof , which are calculable within the iterative 
scheme. The right hand side gives algebraic expressions for the parameters dn and 6„. This method is 
quite powerful and can be used to transform the iterative solution given in terms of Taylor coefficients 
into an equivalent solution in terms of a different set of basis functions than powers of {e.g. those given 
in (jl8p V For a system with non synchronized fiow and initially Gaussian energy distribution, it turns 
out that the first non-vanishing terms already give excellent agreement: 

e-eref wrfiee~''i®-f dae^e"''^®. (19) 

For general values of 77, the expressions are quite cumbersome; at mid-rapidity they reduce to 

di = -^a, (20) 

bi - l + a, (21) 



b 



2a(-15 + 14acr2 + ITcr^a^ + a^) 
81cr2 

298(7^02 - 120a - 165 -I- 166acr2 + iQa'^ + U6a'^a^ 



•3 



6(-15 + Uaa^ + lla^a^ + a^) 



(23) 



with a = a(0o) for ease of the notation. Figure [5] illustrates this approximation. One observes that at 
mid-rapidity a larger value of a(0o) {i.e. a larger initial fiow) leads to negative corrections with respect 
to an initially synchronized flow, meaning a more rapid decrease of central energy density. At large times 
one sees the difference approaching zero, which is a consequence of the energy density approaching zero. 
But if the initial fiow is realistic {i.e. close to synchronized flow, a(9o) ~ 0), the difference decreases even 
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Figure 1: Energy density scaled either by cq = e(6o,?7 = 0) (a) or by eo(To/r)~''/^ = eoexp{ — 19} (b) 
and the difference of flow rapidity to space-time rapidity (c) for a purely longitudinal expansion of an ini- 
tially Gaussian energy distribution. The curves correspond to = ln(r/To) = 0.4, 0.8, 1.2, 1.6, 2.0, 2.4, 3.2 
(in (a) from top to bottom and in (c) counterclockwise) . The curves are the results of the iterative scheme 
up to order 48, and the circles depict a numerical solution of the hydrodynamical equations by a differ- 
encing scheme. 
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faster {i.e. with a larger decay rate in the exponential) than the energy density itself (which decreases 
approximately like the Bjorken energy density cbj = eoe~^^/^). That means that the solution stabilizes 
around the solution with initially synchronized flow, and small perturbations have little effect on the 
long-time behavior. Since for realistic initial conditions (a(0o) close to zero and cr « 4 for LHC) already 
the leading-order correction is very good for all times until freeze-out, one can calculate the time and 
the value of the largest deviation from the reference solution. The maximum of e — e^cf ~ diQe~^'^^ is 
at 9 = l/bi, where the deviation has the value of {di/bi)e~^ . 

A similar analysis can be done to study the dependence of the energy density on the width of the initial 
energy distribution. It is found that the deviation of the energy density e obtained with a Gaussian 
shaped (standard deviation: a) initial energy density from the Bjorken energy density eej can also be 
described with the series ([TS]). In this case, the leading terms are 

esj « dsG^e-"^® + d^O^ e-"-'' . (24) 



6(60,77) 

The normalisation to the initial energy density is necessary to disentangle eff'ects originating from 
gradients from the effect of different initial values at different spatial coordinates. Since the expressions 
for di and 64 at a generic rj are quite extensive this coefficients are only given at rj — 0: 



9( — 2(7^ -I- T]^) 

40-2-45 

'^^^-'^ - -1944^' 

, , 2(-5265 + 428cr2) 

'^^^-'^ - 135(4.2-45) • 

For LHC conditions the leading order term is much more important than the subleading order term 
(except at ry w \/2cr, where it vanishes), i.e. the leading term gives already good agreement. In this case, 
the position and the value of the largest deviation e/e(0o) ~ ercf can be calculated analytically, and one 
finds that the maximum value of (4^2/^2)^^^ taken at 6 = 2/62. In Fig. [31 the approximations are 
plotted at several values of 77. One notes that the next-to-leading order approximation gives excellent 
agreement with the numerical data and the leading-order approximation is also good. 



4.2 Transverse expansion: axisymmetric flow 

Equipped with the finding of the previous subsection we now turn to the axisymmetrical transverse 
expansion. In doing so, we assume for the longitudinal flow 2/(6,77) = 7; consistent with 6(6,77, r) = 
e(6,r) i.e. longitudinal scaling symmetry, which implies an inflnitcly large extension in longitudinal 
direction. The initial transverse energy density distribution is e(6o,7') = 60 exp{— r2/(2K2)} with the 
transverse coordinate r = \/ and k being the standard deviation of the distribution. This transverse 
distribution is certainly to smooth and therefore underestimates the transverse gradients. In this respect 
it may give some lower limit on the transverse flow which is initially zero. Figure [4] shows the evolution of 
the scaled energy density and the transversal flow rapidity for a distribution with initial width k — 7.5 fm 
corresponding roughly to the radius of lead nuclei used for the heavy-ion experiments at the LHC. The 
transversal expansion has a much larger effect on the dynamics than it was the case for the longitudinal 
expansion, since it is not dominated by the initial flow. For instance, at mid-rapidity, the scaled energy 
density drops by ~ 5 % for 6 = ln(T/To) = 2.5 in the case of purely longitudinal expansion(c/. green 
dashed curve in Fig. [SJ, while for the transverse expansion, the scaled energy density is reduced by 58 % 
at the same time (c/. blue dotted curve in Fig. [S]). 

The backbending of the rapidity curve in Fig. 0] (b) for 6 = 2.6 may be attributed to instabilities of the 
truncated iterative scheme for larger times. However, the energy density at the corresponding space-time 
region is already below the freeze-out energy density. 
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Figure 2: Difference of the energy density e of initially not synchronized flow y{Qo) — (l + a(0o))?7 (with 
0(60) = —5 (a) and a(6o) = 5 (b)) to the energy density Ciot of initially synchronized flow (?;(6o) — v)^ 
evaluated at mid-rapidity as functions of the time parameter Q = ln(r/To). In both cases, the initial 
energy density is Gaussian e(9o) = erof(Bo) — eo exp{— 77^/(20-^)} with a = 3.8. The green dashed 
curves are based on (|19p evaluated at leading order, and the blue dotted curves are the corrections 
evaluated at next-to-leading order. The curves are compared to the numerical solution (red solid curves 
with symbols). 
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Figure 3: Difference of the energy density e (normalized to its initial values cq = eo(?7)) of initially 
Gaussian shaped energy density to the Bjorken energy density as a function of = ln(T/To). The 
initial flow is chosen to be synchronized y = rj. The different curves correspond to different longitudinal 
positions (from bottom to top: rj — 0, a, 2a, 3a). The green dashed curves are based on evaluated 
at leading order, and the blue dotted curves are the corrections evaluated at next-to-leading order. The 
red solid curves with symbols depict the numerical solution obtained with a forward integration scheme. 

4.3 Transverse expansion: asymmetric flow 

Transport properties of the sQGP can be extracted from the evolution of initial asymmetries. The most 
intensively studied quantity that encodes such information is the elliptic flow V2, which is used to measure, 
e.g., the shear viscosity of the sQGP. Therefore it is interesting to check whether the iterative scheme 
can cope with asymmetric initial conditions and if the expected change of eccentricity can be observed 
in the results. For this purpose, initial conditions e(6o,X2,a;3) = eo exp{— a;|/(2cr|) — a;§/(2a|)} with 
different widths a2 — 7.5 fm and a^ = 3.75 fm, respectively, along the two transversal axes were used. 
The choice oi a2 ■ a^ — 2 : 1 corresponds to an impact parameter of 6 = 6/5i? w 9fm for lead nuclei (in 
the hard sphere approximation). As in the previous example the initial flow is purely longitudinal with 
y{Qo) = r], and the longitudinal scaling symmetry is kept during the expansion implying independence 
of the energy distribution on the longitudinal coordinate rj. In Fig. [51 the resulting eccentricity 



with e{x2,X3) from a calculation up to order iV = 43 is plotted (solid red curve). 
4.4 Triaxial expansion 

Even triaxial expansion can be treated within the iterative scheme. Since quantities like the elliptic flow 
are measured at mid-rapidity and since in the vicinity of 77 = the energy density does not change much 
in the longitudinal direction, no significant change is expected for the eccentricities at 77 = 0. The dashed 
green curve in Fig. [B] is obtained from a calculation with triaxial initial energy density distribution 
e{Qo, 11^X2, X3) = eoexp{-?72/(2tT2) - a;|/(2o-|) - a:§/(2cr|)} with a = 3.8, a2 = 7.5 fm and (T3 = 3.75 fm 
and initially synchronized Bjorken flow, i.e. 2/(Oo, v) = V ^-nd ^i- = 0- It follows nicely the curve obtained 
for longitudinal scaling expansion, which shows that the spatial eccentricity and thereby the elliptic flow 
at mid-rapidity is not sensitive to the details of the longitudinal flow. However, at > 1.9 instabilities 
due to truncation make a comparison difficult. Curing this either by calculating at higher truncation 
order, optimizing the code or switching to an expansion similar to (llSp remains to be done and is left to 
further studies. 



e = 



/ dx2dx3e{x2,X3){x'^ - xl) 
J dx2dx3e{x2,X3){xl + xD' 



(29) 
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Figure 4: Scaled energy density (a) and transversal flow rapidity ut — arctanh(wT) (b), with vt being 
the velocity in the transversal direction for a transversal expansion of an initially transversal Gaussian 
energy distribution superimposed on the scaling invariant longitudinal flow. The times shown are Q — 
ln(T/To) = 0.4, 0.8, 1.2, 1.6, 2.0, 2.4, 2.6 (from top to bottom in (a) and from bottom to top in (b)). The 
curves are the results of an iterative calculation up to order N = 48. 
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Figure 5: Scaled energy densities at the center xj_ = at mid-rapidity for the Bjorken solution (solid red 
curve) and the longitudinal (green dashed curve) as well as axisymmetric transversal (blue dotted curve) 
and triaxial (orange dash-dotted curve) expansion as a function of time Q — ln(T/ro). The truncation 
orders are N = 48, 48 and 24 for the longitudinal, transversal and triaxial expansion respectively. The 
freeze-out curves correspond to freeze-out temperatures Tf.o. = 165 MeV (left purple dotted line) and 
Ti.o. — 130 MeV (right purple dotted line). The central energy density at initial time is chosen to be 
Co = GSGeVfm^"^, corresponding to LHC energies. 
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Figure 6: Spatial eccentricity at mid-rapidity of the energy density distribution with infinite size in 
longitudinal direction from a calculation up to order A'^ = 43 (red solid line) and the same quantity 
for an additionally longitudinal Gaussian shaped energy density calculated up to order TV = 24 (green 
dashed line). The latter energy density suffers earlier from instabilities and is therefore only depicted 
until e = ln(T/To) = 1.8. 
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Figure 7: Energy density at the edges of the initial distribution. The freeze-out curves and initial 
conditions as well as the truncation orders are the same as in Fig. [5] The longitudinal solution was 
evaluated at ry = cr = 3.8 and the axisymmetric transversal solution at r = k = 7.5 fm. Since the most 
interesting part of the triaxial expansion is the behavior at mid-rapidity it was evaluated at 77 = and 
X2 = and X3 = 3.75 fm {i.e. the steepest point of the initial energy distribution in the transversal plane 
at mid-rapidity). 



5 Limitations 

If it would be possible to calculate the Taylor coefficients up to infinite order, the solution would be 
exact, at least for analytical initial conditions and a smooth equation of state. However, due to limited 
computation power the iteration must be truncated at finite order N. As a result, only a truncated 
series is obtained and errors occur that are at least one order higher than the truncation order. This is 
not necessarily a problem, since in case of heavy-ion collisions one is only interested in a result covering 
a limited time span, for which such corrections may be small. 

In Fig. [5j the energy densities at the center of the fireball for the examples discussed above are plotted 
together with the freeze-out region for typical LHC initial conditions. One observes two facts. First, 
the central energy density of the purely longitudinal expanding system does not differ much from the 
Bjorken solution. This is seen already in Fig. [T](b). Second, the inclusion of transversal expansion leads 
to a much faster cooling of the medium since it expands in more directions. Figure [5] shows that the 
iterative scheme gives reasonable results for the energy density at the center of the fireball. However, 
more problematic than the center are the regions where the gradients are large. To make this statement 
explicit, in Fig. [71 for each of the examples discussed above, the energy density at the edges is plotted. As 
representative for the edge, (one of) the point (s) is chosen, where the initial conditions have the largest 
gradients. One notes that the truncated iterative solutions get unstable quite close to the freeze-out 
region. This is an artefact of the truncated iterative scheme, which demands for calculation with a low 
freeze-out temperature of ^ 130 MeV a higher truncation order. 

It is interesting to find out which properties of the initial condition improve the iterative result and which 
have a negative effect on it. Of course, the result gets better if one can calculate very high orders. To 
achieve this one needs to use initial conditions that can be differentiated most easily. For this reason, 
Gaussian energy distributions were used. An example for initial conditions that are not well suited for 
the iterative method are the ones deduced by Gubser 22 . To understand the computational limits, one 
can estimate the number nk of terms involved in the fcth Taylor coefficient if the initial conditions are of 
Gaussian type. One observes that the Taylor coefficients have the structure of an exponential multiplied 
with a polynomial containing basically all spatial variables and all parameters in all powers up to fc. This 
leads to the estimate, that Uk grows like fc^'i+'^s jyi^ being the number of spatial variables and 7712 
being the number of parameters. The total number of terms involved in the Taylor series truncated at 
order N is therefore ~ Xlf^i fc"i+"2 - j^rm+m^+i ^ 
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Figure 8: Energy density at the center at mid-rapidity for different widths of the initial energy distribu- 
tion. The kink in the curves denotes the time, when the truncated iterative solution gets unstable. The 
initial energy distribution is e(0o) = exp{— ?7^/(2cr^)} with a = 1,2, 4, 8, 16 (from bottom to top). 

A typical operation in the iterative evaluation is the multiplication of two of such series. For such a mul- 
tiplication every term in one series needs to be multiplied with one in the other, leading to ]\j'^i™i+"^2+i) 
terms before simplification. For the transversal asymmetric system discussed in section this results in 
a number of terms which is in the order of 43^ '^ fa 6 x 10^. Therefore, considering additional dimensions 
or letting parameters general leads to a rapidly increasing demand of computation power and reduces 
the calculable order dramatically. A second way to obtain results that are close to the complete solution 
is looking for initial conditions that lead to a quickly converging Taylor series. In this case, less orders 
need to be calculated, which obviously helps the computation. Testing different initial shapes of the 
energy profile we realized that the convergence is better if the initial conditions are more fiat. Steep 
initial profiles lead to rapid oscillations in iterative results which can only be pushed to larger times by 
computing considerably higher orders. A typical example for this behavior is shown in fig. |51 The kink 
that denotes the point, where instabilities of the truncated scheme modify the energy density at the 
central point, occurs the later the broader the initial Gaussian is. 

6 Conclusions 

In this paper we present a scheme that can be used to construct solutions of the equations of relativistic 
fiuid dynamics for differentiable initial conditions. The scheme is based on a formal Taylor expansion of 
the components of the energy-momentum tensor which is used to reformulate the time evolution equations 
as an infinite set of algebraic relations between the Taylor coefficients that can be solved iteratively. In 
the context of AdS / CFT similar methods are used to obtain solutions of the Einstein equations close to 
the boundary of (asymptotic) anti-de Sitter space [371 121] , but to our best knowledge it was not applied 
to hydrodynamics or heavy-ion collisions in the literature so far. The scheme is utilized to address various 
cases of physical interest. It was shown that for the case of purely longitudinal expansion the results 
obtained with the iterative scheme are in excellent agreement with the numerical solution obtained with 
a standard numerical integration routine. One can parametrize the difference between a solution and a 
reference solution in terms of a special set of basis functions {d„0"e~''"®}. The parameters are related 
to Taylor coefficients obtained within the iterative scheme and it was found that already the leading 
order terms give very good agreement with the numerical solution obtained with a standard integration 
scheme. 

For transversal expansion it is possible to find solutions for the axial symmetric as well as for the axial 
asymmetric case. In the latter case, the change of the shape of the energy density due to elliptic flow 
has been considered. Even triaxial expansion needed for more realistic simulations is accessible within 
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the iterative scheme. At mid-rapidity the longitudinal shape of the medium has almost no influence on 
the transversal eccentricity, but at finite rj it is expected to have an effect. 

Concluding we can state that the iterative scheme is a promising tool to investigate the properties of the 
sQGP produced in heavy-ion collision at the LHC via signatures built up during the hydrodynamic part of 
its evolution. With the computational power available for a personal computer at present, its application 
to fluctuating initial conditions or other systems with initially large spatial gradients is somewhat limited. 
Nevertheless there are several applications thinkable, e.g. the calculation of the background field in which 
a jet is moving (c/. the use of hydrodynamics in j34j ) or the emission of real and virtual photons as in 
|35) . Even an application as a benchmark is possible, since the result is exact up to the truncation order 
and should therefore match the results of other hydro-codes. 

A Energy-momentum tensor with derivatives 

Several systems are characterized by an energy-momentum tensor containing derivatives of the basic 
physical degrees of freedom. One important example is viscous relativistic hydrodynamics, e.g. formu- 
lated in the Landau-Lifschitz frame. For such an energy-momentum tensor it is not possible to express 
all components as functions of only. One possible way to cope with that problem is to allow also 
derivatives of T"^ as arguments, i.e. 



with ; denoting the covariant derivative. From (j30|) one can follow the steps outlined in section[2] Another 
possible way would be to construct a similar iterative scheme directly for the basic physical degrees of 
freedom, i.e. in the case of viscous relativistic hydrodynamics for the energy density and the spatial 
components of the four-velocity. 

B Faa di Brunos formula 

The Taylor expansion (I7|) contains time derivatives of Z^", evaluated at the initial time. Since f^''' does 
not depend explicitly on the coordinates the only time dependence originates from the time dependence 
of its arguments T"°. With Faa di Brunos formula [33] it can be decomposed into derivatives of /'"^ 
with respect to T"" and derivatives of T"" with respect to 6: 



To keep the notation as simple as possible, / and g are taken to be scalar functions, but the formula 
can be generalized to g being a vector and / being a tensor function. The set P„ denotes all n-tuples of 
nonnegative integers which fulfill 1 • fci + • • • + n • fc„ — n and D is the differentiation operator. Since m is 
less or equal n, the nth time derivative oi f o g contains only derivatives of g up to order n. Analogously 
the nth time derivative of ff^'^(T°''^) contains only time derivatives of up to order n. 
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